##################
###Reduced Form###
##################

####Baseline
#Subset na omit data
ttt <- na.omit(dis.dat[,c("flus_all","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000")])
#Estimate baseline model
fs1.rf.flu = lm(flus_all ~ comb.np.use_tonnes1000000+t+factor(ccode), data = dis.dat)
fs1.rf.sum.flu <- coeftest(fs1.rf.flu, vcov = vcovCL(fs1.rf.flu, cluster = ~ccode))

####Full
#Subset na omit data
ttt2 <- na.omit(dis.dat[,c("flus_all","flus_all","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])

##Run full model
fs2.rf.flu = lm(flus_all ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = dis.dat)
fs2.rf.sum.flu <- coeftest(fs2.rf.flu, vcov = vcovCL(fs2.rf.flu, cluster = ~ccode))

#####Zero only
##Create a combined meat indicator to assess zero-to-irelevant food production
dis.dat$all_meat_tonnes1000000 <- dis.dat$cattle_meat_tonnes1000000+dis.dat$chicken_meat_tonnes1000000+dis.dat$pigs_meat_tonnes1000000
summary(dis.dat$all_meat_tonnes1000000)
##Subset little-to-irrelevant food production countries
zero.only <- subset(dis.dat, all_meat_tonnes1000000<0.050)
#Remove NAs
zero.only$logGDPpc[is.na(zero.only$logGDPpc)] <- 0
zero.only$comb.np.use_tonnes1000000[is.na(zero.only$comb.np.use_tonnes1000000)] <- 0
#Keep only relevant varialbles 
subset.dat2 <- na.omit(zero.only[,c("flus_all","flus_all","all_meat_tonnes1000000","cattle_meat_tonnes1000000","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","comb.np.use_tonnes1000000","t","lagoutbreak_event", "lagflus_all","logpopulation","logGDPpc","ccode", "year")]) 
#Estimate model
test.iv.rf.flu = lm(flus_all ~ comb.np.use_tonnes1000000+t+lagflus_all+
                  logpopulation+logGDPpc+factor(ccode), data = subset.dat2)
summary(test.iv.rf.flu)
test.iv.rf.sum.flu <- coeftest(test.iv.rf.flu, vcov = vcovCL(test.iv.rf.flu, cluster = ~ccode))
test.iv.rf.sum.flu

#####Rest of the sample
##Subset little-to-irrelevant food production countries
non.zeros <- subset(dis.dat, all_meat_tonnes1000000>0.050)
#Keep only relevant varialbles 
subset.dat3 <- na.omit(non.zeros[,c("flus_all","flus_all","all_meat_tonnes1000000","cattle_meat_tonnes1000000","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","comb.np.use_tonnes1000000","t","lagoutbreak_event", "lagflus_all","logpopulation","logGDPpc","ccode", "year")]) 
#Estimate model
test.iv.rf.2.flu = lm(flus_all ~ comb.np.use_tonnes1000000+t+lagflus_all+
                    logpopulation+logGDPpc+factor(ccode), data = subset.dat3)
summary(test.iv.rf.2.flu)
test.iv.rf.2.sum.flu <- coeftest(test.iv.rf.2.flu, vcov = vcovCL(test.iv.rf.2.flu, cluster = ~ccode))
test.iv.rf.2.sum.flu


#####################################################
#####################################################
###Minimal and Regular/high Production First Stage###
#####################################################
#####################################################

###None production sample
test.iv.di1.flu <- lm(all_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+
                        logpopulation+logGDPpc+factor(ccode), data = subset.dat2 )
summary(test.iv.di1.flu)
test.iv.di.sum1.flu <- coeftest(test.iv.di1.flu, vcov = vcovCL(test.iv.di1.flu, cluster = ~ccode))
test.iv.di.sum1.flu

##Rest of sample
test.iv.di2.flu <- lm(all_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+
                        logpopulation+logGDPpc+factor(ccode), data = subset.dat3)
summary(test.iv.di2.flu)
test.iv.di.sum2.flu <- coeftest(test.iv.di2.flu, vcov = vcovCL(test.iv.di2.flu, cluster = ~ccode))
test.iv.di.sum2.flu


##Export OLS table to Latex
stargazer(fs1.rf.sum.flu,fs2.rf.sum.flu, test.iv.rf.sum.flu, test.iv.rf.2.sum.flu,test.iv.di.sum1.flu,test.iv.di.sum2.flu)
#Obs and diagnostics
stargazer(fs1.rf.flu,fs2.rf.flu, test.iv.rf.flu, test.iv.rf.2.flu, test.iv.di1.flu,test.iv.di2.flu)


#################
#################
###First Stage###
#################
#################

##############
###Baseline###
##############
#Subset NA data
ttt <- na.omit(dis.dat[,c("cattle_meat_tonnes1000000","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000")])
#Beef
fs1.cattle.flu = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+factor(ccode), data = ttt)
#Chicken meat
fs1.chicken.flu = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+factor(ccode), data = ttt)
#Pork
fs1.pigs.flu = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+factor(ccode), data = ttt)

##########
###Full###
##########
##Subset NAs
ttt2 <- na.omit(dis.dat[,c("flus_all","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])

##Beef
fs2.cattle.flu = lm(flus_all ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt2)
##Chicken meat
fs2.chicken.flu = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt2)
##Pork
fs2.pigs.flu = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt2)

##Export OLS table to Latex
stargazer(fs1.cattle.flu, fs1.chicken.flu, fs1.pigs.flu,fs2.cattle.flu, fs2.chicken.flu, fs2.pigs.flu)


